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ABSTRACT 

The 6.7 GHz methanol maser emission, a tracer of forming massive stars, sometimes shows enig- 
matic periodic flux variations over several 10 — 100 days. In this Letter, we propose that this peri- 
odic variations could be explained by the pulsation of massive protostars growing under rapid mass 
accretion with rates of M» > 10~^ Mq yr^^. Our stellar evolution calculations predict that the 
massive protostars have very large radius exceeding 100 Rq at maximum, and we here study the 
pulsational stability of such the bloated protostars by way of the linear stability analysis. We show 
that the protostar becomes pulsationally unstable with various periods of several 10 — 100 days, 
depending on different accretion rates. With the fact that the stellar luminosity when the star is 
pulsationally unstable also depends on the accretion rate, we derive the period-luminosity relation 
log(L/L0) = 4.62 -I- 0.981og(P/100 day), which is testable with future observations. Our models 
further show that the radius and mass of the pulsating massive protostar should also depend on the 
period. It would be possible to infer such protostellar properties and the accretion rate with the 
observed period. Measuring the maser periods enables a direct diagnosis of the structure of accreting 
massive protostars, which are deeply embedded in dense gas and inaccessible with other observations. 
Subject headings: stars: massive — stars: oscillations — masers 



1. INTRODUCTION 

Massive stars have significant impacts on the interstel- 
lar medium through various feedback processes such as 
supernova explosions, stellar winds, and ultraviolet radi- 
ation. These feedback processes would be also important 
for shaping the stellar initial mass function because stars 
arc mostly formed in clusters including massive stars 
(e.g., Lada & Lada 2003). However, our understanding 
of the massive star formation is still limited with obser- 
vational difficulties, one of which arises from the fact that 
forming massive stars are deeply embedded in obscuring 
dense gas. 

Observing maser emission is one of the possible meth- 
ods to see the vicinity (< 10'^ AU) of forming massive 
stars because of the high brightness. It is known that the 
6.7 GHz methanol maser emission can be thought to be 
often associated with circumstellar disks around forming 
massive stars (e.g., Norris et al. 1993; Bartkicwicz ct al. 
2009; Sanna et al. 2010). Some of the methanol masers 
show periodic flux variations over > 10 days (e.g., Goed- 
hart et al. 2004, 2009). Since the 6.7 GHz methanol 
masers are radiatively pumped by infrared emission of 
warm dusts (~ 150 K; Cragg et al. 2005), the periodicity 
could reflect the luminosity variation of nearby forming 
massive stars or accretion disks. 
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Several authors have proposed different explanations 
for this periodic flux variations, e.g., colliding- wind bi- 
nary (van der Walt 2011), and periodic accretion onto bi- 
nary systems (Araya et al. 2010). However, these binary- 
based models do not explain why the periodic variations 
shorter than 10 days have not been found yet, despite 
the fact that a number of OB star eclipsing binaries have 
the periods of < 10 days (Harries et al. 2003; Hilditch et 
al. 2005) 

In this Letter, we propose an alternative picture that 
the periodic variability of the maser emission can be due 
to the pulsation of protostars growing via rapid mass 
accretion with M* > 10"'^ Mq yr~^, which is expected 
in the massive star formation (e.g., Osorio et al. 1999; 
McKee & Tan 2003; Zhang et al. 2005; Behran et al. 
2006; Krumholz et al. 2009). Our previous work shows 
that, by numerically modeling the stellar evolution, the 
massive protostar should have the large radius exceeding 
100 Rq with such the high accretion rate (Hosokawa & 
Omukai 2009; Hosokawa et al. 2010, hereafter HO09 and 
HYOIO respectively). We here examine the pulsational 
stability of the bloated massive protostars with the linear 
stability analysis, and show that the observed periodicity 
of ~ 10 — 100 days can be well explained with the stellar 
pulsation. Our models predict that, with the observed 
period of the pulsation, we could infer the radius, mass, 
and luminosity of the protostars as well as the accretion 
rate onto the protostar. Measuring the maser periods 
would thus make a direct diagnosis of the central small- 
scale (< a few AU) structure of forming massive stars, 
which is beyond observations in the optical and infrared 
bands. 

2. PULSATIONAL INSTABILITY OF ACCRETING 
MASSIVE STARS 

We study the pulsational stability of massive pro- 
tostars growing at constant accretion rates of M* > 
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Fig. 1. — Evolution of the protostcUar radius with various accre- 
tion rates M, = 10"*, 10-=*, 3 X IQ-^, and 6 X 10-=* Mq yr"! 
(taken from HO09). Symbols indicate the stellar models we con- 
duct the linear stability analysis. The circles and stars represent 
the pulsationally stable and unstable models. The shaded layer de- 
notes the instability strip, where the protostar becomes unstable. 

10"'^ Mq yr~^. We here adopt the protostellar mod- 
els with spherical accretion taken from HO09. In Sec- 
tion. 3, we will discuss that the effects of the accretion 
geometry (disk accretion) does not change the qualita- 
tive results from the spherical accretion case. Figure [1] 
shows the evolution of the stellar radius as the stellar 
mass increases with different accretion rates. The out- 
line of the evolution is briefly summarized as follows (see 
HO09 for the details). Initially, the stellar radius grad- 
ually increases with increasing the stellar mass. After 
this stage, e.g., M* > 10 Mq for M* = 10"^ Mq yr-\ 
the protostar contracts by losing its energy via radiation 
(Kelvin-Helmholtz or KH contraction). The stellar cen- 
tral temperature increases during this contraction stage 
and finally reaches lO'' K. The hydrogen burning is ig- 
nited, and the protostar reaches the zero-age main se- 
quence (ZAMS, except with 6 x 10"^ Af© yr^^). After 
this point, e.g., M^ ~ 40 Mq for M» = lO'^ Mq yr'^, 
the stellar radius increases again as the stellar mass in- 
creases. We here note that the maximum stellar radius 
during this evolution is larger with the higher accretion 
rate. This is because the accreting gas has the higher spe- 
cific entropy with the more rapid mass accretion, which 
leads to the higher average entropy in the stellar inte- 
rior (also see IIO09). As a result, the maximum stel- 
lar radius exceeds 100 Rq with the high accretion rates 
M, > 10-3 Mq yr-i. 

We apply the linear stability analysis to the above pro- 
tostellar models (see e.g.. Cox 1980; Unno et al. 1989; 
Inayoshi et al. 2013 for the details). We here consider ra- 
dial (spherical) perturbations with the time dependence 
of e*'^*, where cr = ctr + iai is the eigen frequency, ctr 
is the frequency of the pulsation, and |(Ti| is the growth 
or damping rate of the perturbation. The protostar is 
pulsationally stable (unstable) with the positive (nega- 
tive) ai. If unstable, the perturbation grows until reach- 
ing the non-linear regime, where the pulsation energy is 
dissipated by shock waves near the stellar surface. The 
dissipated energy is partly converted into the kinetic en- 
ergy of periodic outflows (e.g., Appenzeller 1970; Yoon 
& Cantiello 2010). 
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Fig. 2. — The period-luminosity (P-L) relation of forming mas- 
sive protostars. Thin black (blue) curves show the evolution- 
ary tracks with the spherical (disk) accretion at the rates of 
M* = lO—*, 10-3, 3x10-3, and 6x10-3 Mq yr-l. The star sym- 
bols on the tracks denote the pulsationally unstable models. The 
eigen frequencies are plotted for the stable stellar models. The 
shaded layer shows the instability strip as in Figure^] The solid 
line represents the P-L relation given by equation l(y , which fits 
the unstable models. The filled and open symbols represent the 
observed sources, whose distances are measured with the trigono- 
metric parallax and kinematics. The triangles indicate that the 
sources arc associated with ultra/hypercompact Hn regions (Walsh 
et al. 1998, Reid et al. 2009a, Urquhart et al. 2007). 

Symbols in Figure [1] represent the stellar models for 
which we conduct the linear stability analysis. Our cal- 
culations show that the protostar becomes pulsationally 
unstable only when the stellar radius maximally expands 
at a given accretion rate. This instability is caused by 
the K mechanism in the He"*" ionization layer, where the 
radiative energy fiux is blocked and converted into the 
pulsation energy (e.g.. Cox 1980; Unno et al. 1989). In 
the KH contraction stage, the stellar surface tempera- 
ture increases and the He"*" ionization layer disappears. 
The protostar is consequently stabilized and the pulsa- 
tion ceases. Although the protostar would be also unsta- 
ble with the lower accretion rates (< 10"'^ Mq yr-"'^), the 
growth time is much longer than the duration for which 
the star is in the instability strip (cj"^ ^ M^jM^ ~ 10* 
yr) . The perturbation does not grow enough to cause the 
stellar pulsation in this case. The instability strip thus 
docs not extend for M* < 10 M©, where the star be- 
comes unstable with the lower accretion rate (see Figure 

m- 

3. PERIOD-LUMINOSITY RELATION 

Figure [2] shows the evolution of the pulsation period 
and stellar luminosity in the examined cases. In each 
case, the stellar luminosity monotonically increases as 
the stellar mass increases. The pulsation period increases 
with the stellar mass in the early expansion phase, and 
decreases in the later KH contraction phase. The proto- 
star becomes pulsationally unstable around the turning 
point of the period, which corresponds to that of the 
stellar radius as seen in Figure [1] 

Figure [5] also present the observed sources with pe- 
riodic fiux variations in the 6.7 GHz methanol masers, 
whose parameters are summarized in Table [T] The lu- 
minosities of these sources arc estimated with the far- 
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infrared data of Infrared As trono mical Satellite ( IRAS ) 
following ICasoliet alj (|1986D and I Guzman et alj (|2012D . 
We see that the pulsation periods of the unstable models 
are several 10 — 100 days, the same order of the observed 
periods of the maser sources. Although the pulsation pe- 
riod is shorter than 10 days with the low accretion rates 
^* ^ 10""^ Mq yr~^, the instability strip does not cover 
such cases as explained above. This explains well why 
the maser sources which have such the short periodic 
variability have not been observed. 

Our calculations predict that both the period and lu- 
minosity of the pulsationally unstable protostars increase 
with the accretion rate. We thus predict a positive corre- 
lation between the maser periodicity and the stellar lumi- 
nosity. We fit the unstable models by a single power law 
and obtain the period-luminosity (P-L) relation (thick 
hue in Fig. [5]) 

log(i/LQ) = 4.62 + 0.98 log(P/100 days), (1) 

which is similar to that of the Ccpheids (L oc P^^'^), 
which is well- used as a cosmic distance ladder (e.g., Tam- 
mann et al. 2003; Sandage et al. 2004). The above P-L 
relation and the instability strip can explain the distribu- 
tion of the observed parameters in the periodic methanol 
masers within errors of an order of magnitude. We show 
that, in Appendix lA} a similar P-L relation is analytically 
derived. 

The instability strip shown in Figures [T] and [2] is not 
broad. The duration for which a protostar spends in 
the strip, i.e., the prospective lifetime of each periodic 
maser source is only '^ 10^ years. This actually matches 
the rarity expected with observations. The number of 
6.7 GHz methanol masers observed ever is ^ 900 (e.g., 
Pestalozzi et al. 2005; Caswell et al. 2011; Green et 
al. 2012), and 56 masers of that have been monitored 
for more than a year at intervals shorter than a week 
or month C Gocdha r^t al.ll2007l[200llArava et al|[2010l: 
iSzvmczak et al. 20 lit ). About 20 % of such the sources 
show the characteristic periodic variabilities (Table 1). 
Figure [2] might suggest that our models are more ap- 
plicable to the sources with the shorter periods, e.g., 
P < 0.5 year, which is only ~ 5 % of the monitored 
sources. Since the appearance time of the methanol 
masers durin g massive star forrn ation is thought to be 
~ 3 X lO'' yr ()van der Wa"Itll2005| ). the lifetime of the pe- 
riodic maser sources should be as short as expected with 
our calculations. 

An uncertainty of the P-L relation arises from pos- 
sible variations of our stellar evolution tracks. While 
IIO09 study the protostellar evolution with the spheri- 
cal accretion, for instance, HYOIO present that the evo- 
lution slightly changes with different accretion geome- 
tries, e.g., accretion via a geometrically thin disk. We 
also perform the stability analysis of case MD3-D-b0.1 
in HYOIO, whose model has the largest radius among 
the cases for M* = 10"'^ Mq yr""'^. Figure [2] shows that 
the protostar also becomes pulsationally unstable in this 
case, and that the periods are 10 times longer than those 

^ While the "disk accretion" mentioned here assumes the low- 
est entropy of the accreting gas, the rapid mass accretion could 
enhance the entropy (e.g., Hosokawa et al. 2011). In this case, 
the stellar evolution with the spherical accretion would be more 
realistic even if gas accretes through the disk (see HYOIO). 



TABLE 1 

Methanol maser sources with periodic variability 



Source 



Ds, 



Lir 



Ref. 











P 


D 


(Galactic name) 


[day] 


fkpcl 


[10* Lq] 






009.621-1-0.196 


244 


5.2 


30.0 


1 


6 


012.681-0.182 


307 


4.5 


4.2 


2 


7 


012.889-1-0.489 


29.5 


2.34 


2.0 


3 


8 


022.357-f0.066 


179.2 


4.6 


2.4 


4 


9 


037.550-1-0.200 


237 


5.0 


4.5 


5 


9 


188.946-1-0.886 


404 


2.10 


1.2 


1 


10 


196.454-1.677 


668 


5.28 


10.2 


2 


11 


328.237-0.547 


220 


12.0 


117.2 


1 


9 


331.132-0.244 


504 


4.7 


6.2 


1 


7 


338.935-0.062 


133 


2.9 


3.8 


1 


9 


339.622-0.121 


201 


2.6 


1.6 


1 


9 



Note. — Col. 1: name of the source, Col. 2: variability period, 
Col. 3: distance, Col. 4: luminosity estimated with the IRAS data, 
Col. 5 and 6: references of the periods (P) and distances (D). 

Ref erenc es. — (l) [ Goedhart et aT] | |2007I): (2) IGoedhart et al.l 
12004]): O) IGocdhai t et alj 120091): (4 ) ISzvmczak et al.l J2011lh 
(b) lArava et a L 1 2010): ('6') I Sanna et al.l 1200S|): (7) Near kinemati c 
dist anctQ; '(jyiXu ct al.l 1201ll); (9)IGreen fc McClure-Griffithsl 120111) : 
ClO) IReid etai] t2009al ): (iD IHonma et al.l 120071) . 

^Assumed on a fiat rotati on curve with a Galactic circular rotation 
of the Sun, of 246 km s~ IBovv et al.|[2009D . and a solar distance, of 
8.4 kpc (e.g., Rcid et al. 2009b). 



with the spherical accretion at the same rate. This ef- 
fect might explain the observed sources which have the 
longer periods than predicted by the P-L relation ([T]). 

4. CONCLUSION AND DISCUSSION 

In this Letter, we have shown that the pulsation of 
massive protostars could explain the periodic variabil- 
ity of the 6.7 GHz maser sources. Our linear stability 
analysis clearly shows that a rapidly accreting (M* > 
10"'^ Mq yr~^) massive protostar becomes pulsation- 
ally unstable when the star is bloated before reaching 
the ZAMS. Typical periods of the pulsation are sev- 
eral 10 — 100 days, which explain the observed peri- 
odicity well. The period depends on the adopted ac- 
cretion rate, getting longer with the higher rate. On 
the other hand, the protostar with the lower M* be- 
comes unstable but does not produce the periodicity, 
which could also explain why the periodic variability 
shorter than 10 days has not been observed. Our stabil- 
ity analysis predicts the period-luminosity (P-L) relation 
for the pulsationally unstable protostars, log(L/LQ) = 
4.62-f 0.98 log(P/100 day). 

Moreover, our stellar evolution models predict that 
other stellar properties such as the mass, radius, and ac- 
cretion rate as well as the luminosity should depend on 
the pulsation period. We show this in Figure |3l Single 
power-law functions which fit the results are. 
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Fig. 3. — Dependences of the protostellar properties (the stellar 
radius Rt, luminosity L*, mass M*, and accretion rate M*) on the 
period in the stellar-pulsation models. Each line fits the results of 
the stability analysis (symbols) by a single power-law function (see 
equations [l] - |4)l . 

If the P-L relation ([1]) is confirmed by further observa- 
tions, we can infer the above quantities with equations 
([2]) - ([4]) only from the mascr period. As seen in Figure 
[21 the periods of the observed maser sources are longer 
than 10 days. Our stellar-pulsation models predict that, 
in particular to explain the periods longer than 100 days, 
very rapid mass accretion with Af* ^ 3 x 10""^ Mq yr~^ 
is required. In this way, we can make a direct diagnosis of 
the small-scale structure of accreting massive protostars 
and their vicinities (< a few AU), which are difficult to 
see in the optical and infrared bands, by measuring the 
maser periods. 

The stellar pulsation examined above should cause 
the periodic variation of the stellar luminosity X*. Al- 
though it is still uncertain where the maser emission 
is excited around protostars, some observations sug- 
gest that the maser sources are located on circumstcl- 
lar disks at r ~ 10^ AU (e.g., Bartkiewicz et al. 2009). 
The temperature of the circumstellar disk irradiated by 
the stellar radiation obeys Td(r) ~ 120 K (i»/3 x 
10^ L0)i/''(r/lO3 AU)^^/2, where r is the distance from 
the protostar. The variation of the stellar luminosity 
SLf, changes the dust temperature following ST^/T^ ~ 
SLt,/4:L^. Since the response of the dust temperature to 
the irradiating flux variation is so rapid (~ a few sec), 
STd should synchronize with (5i*. As evaluating JL* is 
beyond our linear analysis, we here consider the flux vari- 
ation of Mira variables, which are also excited by the 
K-mechanism. Interestingly, SiO maser sources associ- 
ated with some Mira variables are likely pumped by the 
stellar radiation, and show periodic variations which syn- 
chronize with the stellar pulsation (Pardo et al. 2004). 
With the typical flux variations of ~ 3 magnitudes of 
Mira-type stars (Samus et al. 2009), the amplitude of 
the variation is estimated as ST^/Td ~ 0.3, which re- 
sults in STd ~ 40 K around Td ~ 120 K. Cragg et al. 
(2005) show that the variations of the dust temperature 
6Td > 20 K could produce the observed amplitudes of 
the maser variability {^ 10 — 100 Jy). Therefore, the 
stellar pulsation could produce the flux variations of the 
observed 6.7 GHz methanol masers. 



In the above, we have supposed that the stellar radi- 
ation reaches the irradiated disk surface without signif- 
icant time delays, i.e., the light path from the star to 
the disk surface is optically thin. This should be valid 
with the disk accretion, where the density above the disk 
quickly decreases as the gas falls onto the disk. If out- 
flow is launched from the disk, however, the disk wind 
would enhance the density above the disk and increase 
the optical depth. According to a disk wind model by 
Zhang, Tan & McKee (2013), the optical depth when 
the protostar is pulsationally unstable is estimated as 
~ 5 (M*/10-3 Mq yr-i)0-"(K/3 cm^ g-i), where fi 
is the mean opacity of dusts. The diffusion time with 
this optical depth is ~ 30 days, which is shorter than 
the typical periods of the maser sources. This effect 
would smear out only variabilities shorter than a few 10 
days. Nonetheless, the launching mechanism of the disk 
wind and the detailed density structure above the disk 
are still highly uncertain. Some periodic maser sources 
which have P > 100 days also show rapid fluctuations of 
hght curves over a few 10 days (G009.621, G022.357, and 
G338.935). The above estimate of the diffusion time do 
not explain such the rapid changes. Although the varia- 
tions of these sources might be explained with different 
models rather than ours (also see below), it would be 
also possible that the disk wind is weaker and the diffu- 
sion time is shorter than evaluated above. Simultaneous 
monitoring of the maser sources in the infrared band will 
verify the relation between the periodic variations of the 
stellar luminosity and the maser variability. 

Finally, we compare our pulsation model to alternative 
models for explaining the periodic variability of the 6.7 
GHz methanol maser sources. For instance, van der Walt 
(2011) proposes that radiation from shocked gas in the 
colliding binary wind explains both the periodicity and 
shapes of the light curves of some maser sources. This 
model supposes a particular condition, the binary sys- 
tem surrounded by an ultracompact Hn region (~ 10^ 
AU), and explains the periodic maser sources for which 
ultracompact Hn regions have been observed (triangles 
in Fig [2]). On the other hand, a number of the sources 
presented in Figure [2] do not show associating Hn regions 
at least at 5 and 8 GHz. Our models could explain these 
sources, as the protostar becomes pulsationally unstable 
when the stellar radius is very large and effective tem- 
perature is only Toff ~ 5000 K, with which the UV lumi- 
nosity is too low to create a detectable Hn region. This 
suggests that measuring the UV luminosities of massive 
protostars would be a key for discriminating the models. 
Further observations at the higher frequencies (> a few 
tens GHz) would detect hypercomact Hn regions, which 
are too small and optically thick to be detected at the 
lower frequencies. Such the observations with the higher 
spatial resolution (e.g., with ALMA) are sure to advance 
our understanding of the enigmatic periodic variability 
of the maser sources. 
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APPENDIX 
ANALYTIC DERIVATION OF THE P-L RELATION 

In this Appendix, we analytically derive a P-L relation similar to equation ([T]) with considering the protostcllar 
evolution. Stahler et al. (1986, hereafter SPS86) show that the evolution of an accreting protostar is well understood 
with the balance between following two timescales: the accretion time scale tacc ~ M^/M^, and the Kelvin- Helmholtz 
timescale ^kh = GM^ /L^,R^,. As seen in Figure [TJ the protostar gradually expands in the early stage with iacc < ^kHj 
and turns to contract in the later stage with iacc > ^kh- The turnaround of the radius occurs when iacc — ^kh (also 
see HO09 and HYOIO). According to our linear stability analysis, the protostar becomes unstable at the epoch of the 
maximum radius, i.e., iacc — ^kh (see star symbols in Figure [T|). 

SPS86 also show that the stellar expansion in the early evolutionary stage is well described as i?* oc M,°'^^M^'*^. 
With Kramers' law k cx pT"^/^, which approximates the dependences of opacity in the stellar interior in this early 
stage, the stellar luminosity generally obeys L« a M,^^"_R,~^" (e.g.. Cox & Giuh 1968). Using these relations, we can 

express the luminosity of the unstable protostar (i.e., tacc — ^kh) as a function of the pulsation period P (oc i?* M^, ) 
and obtain the period-luminosity (P-L) relation L, cx P^'^. 
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